#!/usr/bin/python
#coding=utf-8
'''
Analisi finale della figura di interferenza in regime di fotone signolo
---
---
Con un procedimento analogo alla caratterizzazione del piezo vengono calcolate
le funzioni di risposta del piezo utilizzando però i conteggi normalizzati del
PMT.
---
Vengono confrontate le statistiche (gli istogrammi) dei parametri della funzione
di risposta nei due casi (fotodiodo e PMT) usando il test di kolmogorov.
'''
from __future__ import division, print_function
from ROOT import TCanvas, TGraphErrors, TF1, TFile, TTree, TH1D, TLegend
from math import pi
from val2text import tex, set_sig_digit
from myanalisis import set_my_style
import os
from array import array

'''
Definizione dei parametri
'''
#Trasforma la fase del cos^2 nella distanza ottica introdotta 
l = 543 #nm
fase2dist = l/(2*pi)

#Errore sul voltaggio del fotodiodo
svf = 5000/1024

#Grado della polinomiale
n_pol = 2

#Inizializza i file di risultati
path_data = 'data/final/'
path_graph = 'report/graph/final/'
path_table = 'report/table/final/'
path_value = 'report/value/final/'
path_results = 'results/final/'
ofile_name = path_results + 'final.txt'

#importa i coefficienti della funzione di risposta
ifile = open('results/piezo/par_spar.pol2.txt')
pars = []
spars = []
for line in ifile:
    pars.append(float(line.split()[0]))
    spars.append(float(line.split()[1]))

#Inizializza l'albero dei parametri della funzione di risosta a bassa intensità
ofile = TFile(path_results + 'parameters.root', 'recreate')
tree_low = TTree('pars_tree', 'Parametri della funzione di risposta in regime di fotone singlo')
par_pol = array('i', [0])
parameters = array('d', (n_pol + 1)*[0.])
tree_low.Branch('n_pol', par_pol, 'n_pol/I')
tree_low.Branch('pars', parameters, 'pars[n_pol]/D')
par_pol[0] = n_pol + 1

'''
Creazione dei grafici della figura di interferenza senza il limite sul numero
massimo di conteggi 
'''
set_my_style()
for id_file, file_name in enumerate(os.listdir(path_data)):
    #Inserimento dei punti nel grafico
    if not file_name.endswith('dat'):
        continue
    graph = TGraphErrors()
    graph.SetTitle(';V_{piezo} [V];Conteggi normalizzati [ V^{-1} ]')
    with open(path_data + file_name) as file:
        for line in file:
            vf, vp, n = [float(val) for val in line.split()]
            y = n/vf
            sy = ((y/n**.5)**2 + (y/vf*svf)**2)**.5
            if vp > 0.99:
                n = graph.GetN()
                graph.SetPoint(n, vp, y)
                graph.SetPointError(n, 0.00, sy)
    
    #Fit della figura di intrferenza
    fit_func = TF1('fit_func', '[3]*cos(pol2(0))**2 + [4]'.format(4*pi/l), 1.1, 2.1)
    fit_func.SetParameters(0, 3.2, 0.15, 12.5, 0.5)
    graph.Fit(fit_func, 'qm', '', 1.1, 2.1)

    #Grafico
    canvas = TCanvas('canvas', 'canvas', 2)
    graph.Draw('ap')
    fit_func.Draw('same')
    canvas.SaveAs(path_graph + 'figura_di_interferenza_completa_{0}.pdf'.format(id_file))

'''
Creazione dei grafici della figura di interferenza.
'''
for id_file, file_name in enumerate(os.listdir(path_data)):
    if not file_name.endswith('dat'):
        continue
    graph = TGraphErrors()
    graph.SetTitle(';V_{piezo} [V];Conteggi normalizzati [ V^{-1} ]')
    with open(path_data + file_name) as file:
        for line in file:
            vf, vp, n = [float(val) for val in line.split()]
            y = n/vf
            sy = ((y/n**.5)**2 + (y/vf*svf)**2)**.5
            if vp > 0.99 and n < 52000:
                n = graph.GetN()
                graph.SetPoint(n, vp, y)
                graph.SetPointError(n, 0.00, sy)
    
    #Fit della figura di intrferenza
    fit_func = TF1('fit_func', '[3]*cos(pol2(0))**2 + [4]'.format(4*pi/l), 1.1, 2.1)
    fit_func.SetParameters(0, 3.2, 0.15, 12.5, 0.5)
    graph.Fit(fit_func, 'qm', '', 1.1, 2.1)
    for i in range(n_pol + 1):
        parameters[i] = fit_func.GetParameter(i)

    tree_low.Fill()
    prob = fit_func.GetProb()

    #Grafico
    canvas = TCanvas('canvas', 'canvas', 2)
    graph.Draw('ap')
    fit_func.Draw('same')
    canvas.SaveAs(path_graph + 'figura_di_interferenza_{0}.pdf'.format(id_file))

ofile.Write()
ofile.Close()

'''
Compatibilità tra i parametri delle funzioni di risposta ad alta intensità
(ottenute nella caratterizzazione del piezo) e di quelle ad bassa intensità
'''
'''Termine lineare'''
#Parametri degli istogrammi
x_lower_1 = 3.0
x_upper_1 = 3.7
n_bin_1 = 12 
canvas_par1 = TCanvas('canvas_par1', 'canvas_par1', 2)

#Isto termine lineare a bassa intensità
ofile = TFile(path_results + 'parameters.root')
tree_low = ofile.Get('pars_tree')
hist_low_1 = TH1D('pars_low_1', ';Coef. lineare [rad/V];', n_bin_1, x_lower_1, x_upper_1)
tree_low.Project('pars_low_1', 'pars[1]')
hist_low_1.SetLineColor(2)
hist_low_1.SetLineWidth(3)
hist_low_1.SetMarkerStyle(0)
hist_low_1.DrawCopy()

#Isto termine lineare ottenuto dalla caratterizzazione del piezo
file_piezo = TFile('results/piezo/parameters.root')
tree_piezo = file_piezo.Get('par_parerr')
hist_piezo_1 = TH1D('pars_piezo_1', ';Coef. lineare [rad/V];', n_bin_1, x_lower_1, x_upper_1)
tree_piezo.Project('pars_piezo_1', 'pars[1]')
hist_piezo_1.SetLineStyle(2)
hist_piezo_1.SetLineWidth(3)
hist_piezo_1.SetMarkerStyle(0)
hist_piezo_1.DrawCopy('same')

#Legenda
legend = TLegend(0.15,0.92,0.35,0.8);
legend.AddEntry(hist_piezo_1, 'Alta intensita\'')
legend.AddEntry(hist_low_1, 'Bassa intensita\'')
legend.SetFillColor(0)
legend.Draw()
canvas_par1.SaveAs(path_graph + 'par_lineare_isto.pdf')

#Vengono ricalcolati gli istogrammi con numero di bin maggiore per effettuare
#corretamente il test di Kolmogorov
n_bin_1 = 200
#Isto termine lineare a bassa intensità
ofile = TFile(path_results + 'parameters.root')
tree_low = ofile.Get('pars_tree')
hist_low_1 = TH1D('pars_low_1_cum', ';Coef. lineare [rad/V];', n_bin_1, x_lower_1, x_upper_1)
tree_low.Project('pars_low_1_cum', 'pars[1]')
hist_low_1.SetLineColor(2)
hist_low_1.SetLineWidth(3)

#Isto termine lineare ottenuto dalla caratterizzazione del piezo
file_piezo = TFile('results/piezo/parameters.root')
tree_piezo = file_piezo.Get('par_parerr')
hist_piezo_1 = TH1D('pars_piezo_1_cum', ';Coef. lineare [rad/V];', n_bin_1, x_lower_1, x_upper_1)
tree_piezo.Project('pars_piezo_1_cum', 'pars[1]')
hist_piezo_1.SetLineStyle(2)
hist_piezo_1.SetLineWidth(3)

#Compatibilità tra i due istogrammi
prob = hist_low_1.KolmogorovTest(hist_piezo_1)
with open(path_value + 'kolmogorov_prob_par_lineare.txt', 'w') as file_prob:
    file_prob.write(str(set_sig_digit(prob, 2)))

#Calcolo degli istogrammi cumulativi
canvas_par1_cum = TCanvas('canvas_par1_cum', 'canvas_par1_cum', 2)
hist_low_1.ComputeIntegral()
hist_low_1.SetContent(hist_low_1.GetIntegral())
hist_low_1.DrawCopy()
hist_piezo_1.ComputeIntegral()
hist_piezo_1.SetContent(hist_piezo_1.GetIntegral())
hist_piezo_1.DrawCopy('same')
legend.Draw()
canvas_par1_cum.SaveAs(path_graph + 'par_lineare_cum_isto.pdf')

'''Termine Quadratico'''
#Parametri degli istogrammi
x_lower_2 = 0
x_upper_2 = 0.3 
n_bin_2 = 12 
canvas_par2 = TCanvas('canvas_par2', 'canvas_par2', 2)

#Isto termine quadratico a bassa intensità
ofile = TFile(path_results + 'parameters.root')
tree_low = ofile.Get('pars_tree')
hist_low_2 = TH1D('pars_low_2', ';Coef. quadratico [rad/V^{2}];', n_bin_2, x_lower_2, x_upper_2)
tree_low.Project('pars_low_2', 'pars[2]')
hist_low_2.SetLineColor(2)
hist_low_2.SetLineWidth(3)
hist_low_2.SetMarkerStyle(0)
hist_low_2.SetMaximum(4)
hist_low_2.DrawCopy()

#Isto termine quadratico ottenuto dalla caratterizzazione del piezo
file_piezo = TFile('results/piezo/parameters.root')
tree_piezo = file_piezo.Get('par_parerr')
hist_piezo_2 = TH1D('pars_piezo_2', ';Coef. quadratico [rad/V^{2}];', n_bin_2, x_lower_2, x_upper_2)
tree_piezo.Project('pars_piezo_2', 'pars[2]')
hist_piezo_2.SetLineStyle(2)
hist_piezo_2.SetLineWidth(3)
hist_piezo_2.SetMarkerStyle(0)
hist_piezo_2.DrawCopy('same')

#Legenda
legend = TLegend(0.15,0.92,0.35,0.8);
legend.AddEntry(hist_piezo_2, 'Alta intensita\'')
legend.AddEntry(hist_low_2, 'Bassa intensita\'')
legend.SetFillColor(0)
legend.Draw()
canvas_par2.SaveAs(path_graph + 'par_quadratico_isto.pdf')

#Vengono ricalcolati gli istogrammi con numero di bin maggiore per effettuare
#corretamente il test di Kolmogorov
n_bin_2 = 200
#Isto termine quadratico a bassa intensità
ofile = TFile(path_results + 'parameters.root')
tree_low = ofile.Get('pars_tree')
hist_low_2 = TH1D('pars_low_2_cum', ';Coef. quadratico [rad/V^{2}];', n_bin_2, x_lower_2, x_upper_2)
tree_low.Project('pars_low_2_cum', 'pars[2]')
hist_low_2.SetLineColor(2)
hist_low_2.SetLineWidth(3)

#Isto termine quadratico ottenuto dalla caratterizzazione del piezo
file_piezo = TFile('results/piezo/parameters.root')
tree_piezo = file_piezo.Get('par_parerr')
hist_piezo_2 = TH1D('pars_piezo_2_cum', ';Coef. quadratico [rad/V];', n_bin_2, x_lower_2, x_upper_2)
tree_piezo.Project('pars_piezo_2_cum', 'pars[2]')
hist_piezo_2.SetLineStyle(2)
hist_piezo_2.SetLineWidth(3)

#Compatibilità tra i due istogrammi
prob = hist_low_2.KolmogorovTest(hist_piezo_2)
with open(path_value + 'kolmogorov_prob_par_quadratico.txt', 'w') as file_prob:
    file_prob.write(str(set_sig_digit(prob, 2)))

#Calcolo degli istogrammi cumulativi
canvas_par2_cum = TCanvas('canvas_par2_cum', 'canvas_par2_cum', 2)
hist_low_2.ComputeIntegral()
hist_low_2.SetContent(hist_low_2.GetIntegral())
hist_low_2.DrawCopy()
hist_piezo_2.ComputeIntegral()
hist_piezo_2.SetContent(hist_piezo_2.GetIntegral())
hist_piezo_2.DrawCopy('same')
legend.Draw()
canvas_par2_cum.SaveAs(path_graph + 'par_quadratico_cum_isto.pdf')

ofile.Close()
file_piezo.Close()
